Convergent Cross Mapping (CCM)
School on Physics Applications in Biology, January 2018
This tutorial starts by presenting the general idea of Convergent Cross Mapping (Ye et al. 2017), next it shows some applied examples using syntetic data. In the last section a real data analysis is proposed as an exercise. For these activities, you will need the most recent version of R and the rEDM package installed in your working computer.
Convergent Cross Mapping: From Chaos to Chaos-ality
Syntetic Data
First, We will need to simulate a deterministic discrete-time dynamics with chaotic behavior. To create a time series from this dynamics with 150 time steps run the commands below in the R console:
## Two vectors to store data
X <- c()
Y <- c()
## Initial values
X[1] <- 0.1
Y[1] <- 0.1
## Iterate the dynamics 150 time steps
for(i in 2:150){
X[i] <- 3.77*X[i-1]*(1-X[i-1])
Y[i] <- 3.82*Y[i-1]*(1-Y[i-1]-0.05*X[i-1])
}
XY<-as.data.frame(cbind(X,Y))Here we have the time series of two coupled species in a chaotic regime. In this first example the variable X is causing changes in Y. However, populations do not seem to be related.
The correlation coeficient between X and Y is shown in the figure below.
Although the data comes from a model with causality, X and Y are not correlated. This is a clear example that “Lack of correlation does not imply lack of causation”.
Cross Mapping1
How can we extract this causality from the coupled dynamics? As we have seen in a previous section, a generic property of lagged-coordinate embedding is that points x(t) on MX map 1:1 to points m(t) on M and local neighborhoods on Mx map to local neighborhoods on M.
It follows that for two variables X and Y that are dynamically coupled, local neighborhoods on their respective lagged reconstructions, MX and MY, will map to each other since X and Y are essentially alternative observations of the common attractor manifold M.
Convergent cross mapping determines how well local neighborhoods on Mx correspond to local neighborhoods on My and vice versa. To do so, a manifold Mx is constructed from lags of variable X and used to estimate contemporaneous values of Y. Similarly, a manifold My is constructed from lags of variable Y and used to estimate contemporaneous values of X.
To do so, we first need to obtain the optimal embedding dimension for both variables using the simplex function.
options(warn = -1)
simplex_X<-simplex(X,silent=T)
plot(simplex_X$rho,type='o')
E_star_X<-which.max(simplex_X$rho)
print(paste('E*(X) =',E_star_X))## [1] "E*(X) = 2"
simplex_Y<-simplex(Y,silent=T)
plot(simplex_Y$rho,type='o')
E_star_Y<-which.max(simplex_Y$rho)
print(paste('E*(Y) =',E_star_Y))## [1] "E*(Y) = 2"
Cross-mapping from Mx to My (X_xmap_Y)
Now that we have the optimal embeddings, we can construct the shadow manifolds Mx and My using the make_block function from rEDM package. Use this command to download the function code and load it in your R workspace (it will be included in rEDM itself in future releases):
Now for the shadow manifold:
Shadow_MXY<-make_block(XY,max_lag = 2)
Shadow_MX<-Shadow_MXY[,2:3]
Shadow_MY<-Shadow_MXY[,4:5]
head(Shadow_MXY)## time X X_1 Y Y_1
## 1 1 0.1000000 NA 0.1000000 NA
## 2 2 0.3393000 0.1000000 0.3418900 0.1000000
## 3 3 0.8451417 0.3393000 0.8373481 0.3418900
## 4 4 0.4934071 0.8451417 0.3851034 0.8373481
## 5 5 0.9423361 0.4934071 0.8682788 0.3851034
## 6 6 0.2048571 0.9423361 0.2806179 0.8682788
To better understand the cross-mapping, lets start by performing one single prediction using the index of a single predictor.
Here, we are using the therm “prediction”, but instead of predicting future values of X, we will “predict values of Y(t) using X(t) and vice versa. This”cross-prediction" is performed between different variables but for the same point in time.
Next, we obtain the indexes of the E*+1 nearest neighbors from the given predictor generating the simplex_Mx.
dist.matrix_X <- as.matrix(dist(Shadow_MX, upper=TRUE))
neigb_X <- order(dist.matrix_X[predictor,])[2:4]
neigh_X_print<-c(neigb_X)
print(paste('simplex_Mx:',list(neigh_X_print)))## [1] "simplex_Mx: c(31, 93, 130)"
The following plot presents Mx with the predictor (red dot) and respective simplex_Mx (blue dots).
The cross-mapping process starts by finding (mapping) the simplex_Mx in My, creating the simplex_My. Note that simplex_My has the same indexes as simplex_Mx. The figure below shows My and simplex_My (green dots).
The simplex(My) will then be used to estimate the contemporaneous value of the predictor in My, obtaining the predicted value Y(tilda).lib<-lib <- c(1, NROW(Shadow_MXY))
block_lnlp_output_XY <- block_lnlp(Shadow_MXY, lib = lib, pred = lib, columns = c("X",
"X_1"), target_column = "Y", stats_only = FALSE, first_column_time = TRUE)
observed_all_Y <- block_lnlp_output_XY$model_output[[1]]$obs
predicted_all_Y <- block_lnlp_output_XY$model_output[[1]]$pred
pred_obs_Y<-as.data.frame(cbind(predicted_all_Y,observed_all_Y))
colnames(pred_obs_Y)<-c('Predicted Y','Observed Y')
head(pred_obs_Y)## Predicted Y Observed Y
## 1 0.4975515 0.8373481
## 2 0.6043905 0.3851034
## 3 0.7493451 0.8682788
## 4 0.5365754 0.2806179
## 5 0.6601371 0.7601691
## 6 0.5649814 0.6072697
Taking every prediction in the series will provides us with many comparisons between predicted and observed. This is better vizualized in the plot below. The red dot represents the predictor used in the example above.
Cross-mapping from My to Mx
Similarly to the previous mapping, we obtain the indexes of the E*+1 nearest neighbors from the given predictor. However, since we are cross-mapping from My to Mx, the simplex_My is generate in shadow maninfold My.
dist.matrix_Y <- as.matrix(dist(Shadow_MY, upper=TRUE))
neigb_Y <- order(dist.matrix_Y[predictor,])[2:4]
neigh_Y_print<-c(neigb_Y)
print(paste('simplex_My:',list(neigh_Y_print)))## [1] "simplex_My: c(135, 5, 54)"
The following plot presents My with the predictor (red dot) and respective simplex_My (blue dots).
Next, we map the simplex_My in Mx, creating the simplex_Mx. Analogously, simplex_Mx has the same indexes as simplex_My. The figure below shows Mx and simplex_Mx (green dots).
The simplex(Mx) will then be used to estimate the contemporaneous value of the predictor in Mx, obtaining the predicted value X(tilda).
lib<-lib <- c(1, NROW(Shadow_MXY))
block_lnlp_output_YX <- block_lnlp(Shadow_MXY, lib = lib, pred = lib, columns = c("Y",
"Y_1"), target_column = "X", stats_only = FALSE, first_column_time = TRUE)
observed_all_X <- block_lnlp_output_YX$model_output[[1]]$obs
predicted_all_X <- block_lnlp_output_YX$model_output[[1]]$pred
pred_obs_X<-as.data.frame(cbind(predicted_all_X,observed_all_X))
colnames(pred_obs_X)<-c('Predicted X','Observed X')
head(pred_obs_X)## Predicted X Observed X
## 1 0.7513307 0.8451417
## 2 0.4587703 0.4934071
## 3 0.9306259 0.9423361
## 4 0.2179312 0.2048571
## 5 0.6259113 0.6140977
## 6 0.9040863 0.8934210
Below is the plot for all predictions and observations as well as the Pearson’s correlation coefficient. The red dot represents the predictor used in the example above.
What does it all mean?
From the previous results, if we use Y to estimate X we obtain good predictions. Since there is information about X inside Y, we can use Y to estimate X. In another words, if X causes Y, the cross mapping skill from My to Mx will generally gives us good results (high correlations). Important to note here is the inverse relation between cross mapping and causality.
Convergent part of CCM
Convergence means that cross-mapped estimates improve in estimation skill with time-series length L (sample size used to construct a library).
With more data, the trajectories defining the attractor fill in, resulting in closer nearest neighbors and declining estimation error (a higher correlation coefficient) as L increases. Thus, CCM becomes a necessary condition for causation.
Exercises
What happens when there is a invertion in cause-effect realation? Change the model so the Y variable is causing changes in X.
What if there is no interaction between X and Y?
What about both variables interacting to each other?
Identifying causal networks is important for effective policy and management recommendations on climate, epidemiology, financial regulation, and much else. In the following exercise, you should use CCM to identify causality between anchovy landings in California and Newport Pier sea-surface temperature.
Learn more
- Sugihara G. and R. M. May. 1990. Nonlinear forecasting as a way of distinguishing chaos from measurement error in time series. Nature 344:734–741.
- Sugihara G. 1994. Nonlinear forecasting for the classification of natural time series. Philosophical Transactions: Physical Sciences and Engineering 348:477–495.
- Anderson C. & Sugihara G. Simplex projection - Order out the chaos. http://deepeco.ucsd.edu/simplex/
- “Simplex projection walkthrough”, a tutorial by Owen Petchey.
This section is adapted from: Sugihara et al., Detecting Causality in Complex Ecosystems (Supplementary Materials), Science, 2012.↩